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SUMMARY 


This  fifth  and  final  volume  of  the  report,  "Hydrodynamic  Effects  of 
Nuclear.  Explosions,  "  presents  new  theoretical  developments  for  two 
problems.  The  first  is  the  determination  of  the  waves  resulting  from 
the  passage  of  a  high-pressure  disturbance  over  the  free  surface  of  a 
body  of  water.  This  would  occur  in  the  case  of  a  burst  on  or  over  land 
near  a  shore  and  is,  therefore,  of  interest  to  the  Five  City  Study  in 
which  three  nuclear  surface  bursts  are  near  rivers  or  bays. 

The  second  topic  is  motion  of  the  ground  water  table  induced  by  a  surface 
burst.  This  problem  is  of  interest  for  the  determination  of  the  migration 
of  radioactive  contaminants  and  is  also  applicable  to  three  cities  of  the 
Five  City  Study. 

Only  theoretical  development  is  given  here.  Typical  methods  of 
application  may  be  found  in  Volume  IV  of  this  report  subtitled  "Five 
City  Study .  " 
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NOMENCLATURE 
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F(a,  b,  c,  z) 
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I 
J 
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L 

P 

P 

o 

P 

n 

P* 

IP 

Q 

Q 


decay  constant  in  error  function 
Laplace  transform  function 
friction  factors 
acceleration  of  gravity 
water  depth 

variable  in  Hankel  transform 
summation  indices 
variable  in  Laplace  transform 
time 

horizontal  particle  velocity 
horizontal  distance  from  shore 
error  function 

complimentary  error  function 
arbitrary  constant 

Laplace-Hankel  transform  of  dimensionless  surface 
elevation 

hypergeometric  function 
function 

Heaviside  step  function 
integral 

Bessel  function  of  nth  order 
length  scale 

pressure  acting  upon  surface 
Magnitude  of  pressure  step 
Legendre  function 
pressure  front 

Laplace  transform  of  dimensionless  pressure 
Laplace  -  Hankel  transform  dimensionless  pressure 
magnitude  of  particle  velocity  vector 
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Part  II 


d 

g 

h 


variable 
time  scale 

horizontal  velocity  of  pressure  front 
wind  velocity 
beach  slope 

exponential  decay  factor 
dimensionless  pressure 
Dirac  delta  function 
dimensionless  surface  elevation 
dynamic  response 
surface  elevation 
T 

fluid  density 
air  density 
dimensionless  time 
shear  stress  at  surface 
shear  stress  at  bottom 
transcendental  function 
contour  of  integration 
gamma  function 
variable  in  transform 
refer  to  partial  derivatives 
refers  to  peak  depression 
refers  to  peak  elevation 


distance  between  well  and  reference  point 
various  functions  as  required  in  analyses 
water  depth 
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NOMENCLATURE 

(Concluded) 


c 

e 

x 

V 

p 

p 

<p 

7 

v2 

0* 

«o 

K 

0n 

etc. 


dimensionless  water  level  ratio 
angle 

distance  between  crater  and  river 

variable 

fluid  density 

H/R 

potential  function 

gradient  operator 

Laplacian  operator 

refers  to  dimensionless  variables 

refers  to  initial  value 

refers  to  transient  water  level 

refers  to  stream  water  level 

refers  to  successive  terms  in  series  solution 

refers  to  partial  derivatives 


NOMENCLATURE 

(Continued) 


k  =  constant  in  Darcy  law 
m  =  porosity  of  soil 
n  =  summation  index 
p  =  pressure 
q  =  unit  discharge 
r,  0,  z  =  cylindrical  coordinates 
s  =  drawdown  of  well 
t  =  time 

u  =  particle  velocity  component  in  r-direction 

w  =  particle  velocity  component  in  z-direction 

w,  x,  y,  z  =  dummy  variables  as  required 

Ap  =  particle  radius  of  bed  material 

A  ,  B  =  constants 
n  n 

C  =  concentration 
C,  D  =  constants 
D  =  Domain 

Dx,  Dy  =  dispersion  coefficients 
E,  F  =  variable  functions 
H  =  height  of  fluid  at  r  - 

Iq,  Kq  =  Bessel  functions  of  1st  and  2nd  kinds,  respectively 
K  =  permeability 
Q  =  discharge  rate 
R  =  cavity,  crater  or  well  radius 
Sj  =  free  surface 
T  =  hydraulic  transmissivity 
U  =  velocity 
a,  0  =  variables 

y  =  specific  weight  of  fluid 
V  =  datum 

C  =  specific  yield  =  effective  porosity 


PART  I 

WAVES  GENERATED  BY  A  TRAVELING  DISTURBANCE 


1.  ESTABLISHMENT  OF  THE  BASIC  EQUATIONS 


Consider  a  two-dimensional  pressure  disturbance  arriving  at  the  free 
surface  of  a  body  of  water,  as  illustrated  in  Fig.  1.  The  response  of 
the  water  to  such  a  disturbance  will  be  analyzed  by  means  of  the  long 
wave  equations  valid  for  shallow  water: 

Continuity: 


n,  + 


[u  (h 


r,)] 


=  0 


Momentum: 


(1) 


u.  +  uu 

t  X 
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p 
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(2) 


where 

x 

t 

h 

V 
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g 

P 

P 


distance  from  shore 

time  after  arrival  of  pressure  at  shore 
h(x)  =  still  water  depth 

surface  elevation  around  still  water  level 

horizontal  component  of  water  particle  velocity,  assumed 
constant  over  a  vertical 

acceleration  due  to  gravity 

water  density 

pressure  acting  on  the  surface 
shearing  stress  at  the  surface 
shearing  stress  at  the  bottom 


and  subscripts  are  used  to  denote  partial  differentiation  with  respect  to 
themselves. 


Of  couse,  it  is  possible  to  treat  these  equations  successfully  by  numerical 
techniques  (finite  differences,  the  method  of  characteristics).  However, 
through  suitable  approximation,  it  is  possible  to  derive  analytic  solutions, 
in  certain  instances,  which  are  sufficient. 


p(*,t) 


Figure  1 

Problem  configuration 
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The  shearing  stresses,  t  and  t,  ,  may  be  written 

S  D 

%  ■  Oa  UI 

Tb  =  Pf2u2 

where  p^  is  the  air  density  and  is  the  wind  velocity  near  the  surface. 
Although  may  be  large,  and  are  small  so  that  Tg  may  be 

neglected.  Similarly,  although  p  is  large,  f^  and  u  are  small  so  that 
may  be  neglected.  Hence,  the  shearing  stress  term  will  be  neglected; 
a  good  approximation  except  near  the  shore  where  h  +  Tj  =  0, 

The  equations  may  be  linearized  by  assuming  the  convective  inertia  term, 
uu^,  to  be  small  and  that  t-  may  be  neglected  in  comparison  with  h  (a 
crude  approximation  near  the  shore).  Then  the  equations  become: 

Continuity: 


\  +  (uh)x  =  0 


(3) 


Momentum: 


U  =  -  g  7] - P 

t  ”  *X  P  X 


Differentiating  Eq.  3  by  5/Bt  and  Eq.  4  by  ?/Bx  gives 


TJ. .  +  u  h  +  u  h  =  0 
tt  t  x  xt 


and 
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xt 


gn 


XX 
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p 

XX 


(4) 


(5) 


(6) 


3 


Eliminating  and  u^  from  Eq.  5  by  means  of  Eq.  4  and  Eq.  6  results 
in  an  equation  for  Tj: 


7]..  -  h  (grj  +  4  p  )  -  h  (grj  +  4  P  )  =  0 
't t  x  B  'x  p  x  6  'xx  p  xx 


(7) 


For  a  uniformly  sloping  beach,  h  =  ax,  one  has 


XT)  +  7]  -  —  T). . 
'xx  'x  ag  'tt 


-  ~  (xP  +  P  ) 

Pg  XX  X 


(8) 


By  defining  time  and  length  scales 


T  =  U/ag 


L  =  U^/ag 


(9) 

(10) 


where  U  is  some  (constant)  typical  velocity,  one  may  introduce  the 
dimensionless  variables 

C  *  ri/L 

\  =  x/L 
T  =  t/T 
y  2  p/pgL 


(ii) 


j 


which,  when  substituted  in  Eq.  8,  gives 


y  r  +  £  -  C  =  -  (x  y  +  y  ) 
*  XX  X  tt  '  \x  'x 


(12) 


Assuming  no  initial  disturbance,  the  initial  conditions  are 


C  (\.0)  =  CT  (\.0)  =  o 


(13) 
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The  boundary  condition  is  imposed  that 

C  -  0  as  \  -  ®  (14) 

The  general  solution  of  this  initial  boundary  value  problem  can  be  found 
by  first  applying  the  Laplace  transform  to  t,  making  the  change  of 
variable 


?  =  yv 

and  then  applying  the  Hankel  transform  of  zeroth  order  to 
The  Laplace  transform  of  Eq.  12  with  respect  to  T  is 
\  f"  +  f'  -  s2f  =  -  (\  F"  +  F') 
where 

X 

f(\.  s)  =  ^  C(x»  t)  e’ST  dr 

ff(\,  s)  =  \  y(x,  t)  e‘sT  dT 

'0 

and  primes  denote  differentiation  with  respect  to 
Introducing  the  change  of  variable 

-•  ■  VT 

Eq.  15  becomes 

d2f  1  df  .  2f  fdZW  ,  1  dP 

1?  +  t^‘4s  f  =  '^  +  Td^ 


(15) 


(16) 

(17) 


(18) 


(19) 


5 


Letting 


< 


F(k)  = 

r®  • 

\  f(S)  t  J  (k»)  d» 

J0  0 

Q(k)  = 

00 

[  P(»)  '  J  (kf)  d: 

Ja  0 

and  using  the  property 


-  k2  F 


the  Hankel  transform  of  Eq.  19  gives 


F 


k2_0_ 
k2  +  4s2 


(20) 

(21) 


(22) 


(23) 


Then  the  solution  obtained  by  inversion  transforms  is  found  to  be 

c+i°°  ■»  2 

C(x.r)  =  -V(X.T)  +  25i-  \  ,  d.  eST  \  dk  k  J0(k-)  J*  ° 

'c-i'”  0  (k  +  4s  ) 


.(24a) 


The  first  part  of  this  expression  clearly  represents  the  hydrostatic 
response,  whereas  the  integral  term  represents  the  dynamic  response 
which  will  be  denoted  by  That  is, 


C(\»  r)  =  -  y(\,  r)  +  Cd 


(24b) 


Before  presenting  detailed  analyses  of  this  solution  in  specific  cases,  it 
is  of  interest  to  present  a  simplified  calculation  of  the  small  time  response. 
Neglecting  the  gravity  terms  in  Eq.  7  results  in 


e 


TT 


X  7 XX  +  y 


X 


(25) 
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A  pressure  wave  possessing  a  sharp  front  moving  at  speed  U  in  the 
positive  \  direction  can  be  represented  as: 

y(\,  r)  =  y(r  -  \)  =  h(t  -  \)  P*  (r  -  \)  (26) 


where  H  is  the  Heaviside  step  function  and  P*  is  continuous.  Because 
of  the  form  of  P*,  one  may  replace  dy/^\  by  -  (dy/br)  in  Eq.  25: 


C  %  v  y  -  v 

STT  N  TT  T 


Making  use  of  Eq.  13  and  the  identity 

T  T~\ 

[  H(t  -  \)  P:':  (t  -  \)  dT  =  H(r  -  \)  \  P-(T)d; 

'o  '0 

one  obtains 


(27) 


(28) 


C 


f  pT-\  1 

H(r  -  X)  \  P*  (T  *  \)  -  \  P ••••(?)  d* 

L  .>0  J 


r  «  1 


(29) 


This  may  be  generalized  for  an  arbitrary  bottom  y  =  -h(\)  and  any 
pressure  distribution  traveling  unchanged  at  constant  speed,  y  -  y(t  -  \): 


rT 

C  *=  h(\ )  y(r  -  \)  -  h  (\)  •>  (r  -  \)  dr  r  «  1  (30) 

*  ''O 


If  v  has  a  sharp  front  ahead  of  which  the  pressure  is  zero 

,T-\ 


H(t  -  \)  {h(\ )  P:;:(T  -  \)  -  h^(X)  ^  P*(»)  d?} 


T  «  1 


(31) 


These  expressions  are  a  first  approximation  only  but  are  useful  since  P* 
and  h  may  be  quite  arbitrary.  Improvements  of  a  higher  order  can  be 
obtained  by  investigating  the  full  boundary  value  problem. 
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2.  STEP- FUNCTION  PRESSURE  OVER  A  UNIFORMLY  SLOPING  BEACH 

The  simplest  model  of  the  pressure  assumes  a  sharp  front  moving  at 
unit  velocity,  the  pressure  being  zero  ahead  and  unity  behind.  In  other 
words 

y(\,  T )  =  H(r  -  \)  =  H(t  -  ?2)  (32) 

where  H  is  the  Heaviside  step  function.  The  sharp  front  of  this  model 
is  a  realistic  representation  of  the  traveling  shock,  although  the  pressure 
behind  the  shock  will  actually  decay.  Such  a  decaying  pressure  will  be 
treated  in  Section  3. 


The  composite  Laplace- Hankel  transform  of  y  is 


,  1  /  k2' 

Q(k,sl  =  —j  exp  v- 

2s 


(33) 


Inserting  this  expression  into  Eq.  24  it  is  found  that 


I  k2^ 

r*  1  rc4l0r  exP  vST  -  4 sJ 

r.  .  =  2  \  k  J  (k •)  —T-  \  — U - 

d  -  °  2771  Jc-i-  k2  +  4s2 


ds  dk 


J0 


(34) 


We  may  evaluate  the  Laplace  inversion  integral  first  and  then  deal  with 
the  Hankel  integral.  Although  the  behavior  of  the  Laplace  integral  for 
large  t  can  be  found  by  the  method  of  steepest  descent,  as  will  be  dis¬ 
cussed  in  Section  5,  the  results  are  not  too  significant,  physically.  Here 
we  shall  be  concerned  with  the  small  time  response  only. 


The  response  at  small  t  is  represented  by  the  transform  at  large  s. 
Therefore,  we  expand  the  integrand  for  large  s: 


I  % 


1 

Ziri 


,c+i<» 

c-i® 


l 


(-k2) 


n-  1 


i  a  2.n 
(4s  ) 


exp 


(35) 
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Each  integral  in  the  expansion  may  be  evaluated: 


iE 

c  n*l 


(-1) 


n- 1  r 


n-  1/2 


2n-  1 


(k  V75 


(36) 


so  that: 


-£ 


(_un-l  Tn-l/2 


S”Jo,k’>  J2„-l(k  T/Tl  dk 


(37) 


The  integral  factor  here  may  be  evaluated  in  terms  of  Gauss'  hyper¬ 
geometric  function: 

J  J0(k5)  J2n_j  (ky?)  dk  =  r‘1/2  F  (n;-n+l,  1,*)  X<T 

=  0  \  >  T 

which,  in  turn,  may  be  expressed  in  terms  of  Legendre  polynomials: 

LJo(k?)  J2n-1  (kV^>  dk  =  H  |T  -  X)  Pn.!  (l  -  2*)  (38) 

Inserting  this  result  into  Eq.  37  gives 


=  0 


\  <T 


X  >  T 


(39) 


Combining  this  with  the  static  component  -y  =  -  H  (r  -  \)  gives 


c .  { y  (-T>np  „(*  -2?)}  h<t- 

n^  1 


X) 


T  «  1 


(40) 
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This  series  may  be  summed  explicitly,  giving: 

c  %  (l  +  2r  +  r2  -  4xf 1/2  -  l  \  <  r 

(41) 

'=0  X  >  T 

Since  r  and  \  are  very  small  compared  to  unity,  this  expression  may  be 
approximated  as 


C  %  2y  -  t 

=  0 


which  is  seen  to  be  the  first  term  of  the  series  in  Eq.  40,  the  succeeding 
terms  being  small.  As  illustrated  in  Fig.  2,  the  free  surface  is  deformed 
only  beneath  the  pressure,  having  a  sharp  front  at  the  pressure  front. 

The  mathematical  profile  OA'CD  satisfies  conservation  of  mass.  Due  to 
the  presence  of  the  bottom,  however,  the  physical  system  cannot;  perhaps 
the  new  shoreline  at  A  should  be  interpreted  as  a  sink-like  singular 
point.  Singular  behavior  at  a  shoreline  has  been  investigated  by  Ho  and 
Meyer  (1962). 


Reverting  to  physical  variables,  it  is  found  that  the  peak  elevation,  r^, 
and  the  peak  depression,  rjp,  are 

Pat 

r  _ 

'E  ~  pU 

P  rv2  t 
o 

''D  *  "  pU(2  +  a) 

where  P^  is  the  magnitude  of  the  pressure  step. 

This  case  may  be  extended  to  consider  a  step  of  finite  duration.  The 
pressure  may  be  thought  of  as  a  superposition  of  two  step  functions: 

y  (\,T)  =  H(t  -  X)  -  H(r  -  Tq  -  \)  (44) 


11 


Analysis  similar  to  the  foregoing  yields  the  result  illustrated  in  Fig.  3; 
the  length  AE  depends  on  T  and  Tq  and  is  an  almost  uniform  depression. 
Again,  the  maximum  water  elevation  is  =  T. 
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Figure  3 
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Finite  duration  step-pressure  response 
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3.  DECAYING  PRESSURE  OVER  A  UNIFORMLY  SLOPING  BEACH 


In  this  second  model,  we  assume  a  pressure  distribution  that  decays 
behind  a  sharp  front,  as: 


y{\,  T )  =  H  (r  -  Y)  e"^T'^ 


(45) 


Then  we  have 

Q  (k,s)  = 


2s  (s  +  0) 


(46) 


The  first  part  of  this  expression  is  the  same  as  that  of  Eq.  33  and  so 

represents  the  dynamic  response  of  that  model,  which  we  may  denote 

.  „(1) 

by  ,  enabling  us  to  write 


.01 


-  0 


f 

Jr 


T':) 


dr* 


by  the  convolution  theorem. 


(47) 


Again  concentrating  on  behavior  at  initial  moments  (t  «  1)  we  may  use 
the  small-time  approximation  of  i.  e.  ,  Eq.  39  to  evaluate  the  con¬ 

volution  integral: 


(\»  T*) 


-0  (t-t::) 


d  T* 


CT 

=  -  0H  (r  -  Y)  \  (l  +  It"  +  r 
Y 


4\)'1/2  e 


0(t-t*) 


dr* 
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i 
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Expanding  the  radical  and  neglecting  terms  above  first  order  (since 
T  «  1 )  gives 


+  2\  -  T  -  \  e 


■/3(t-\) 


0(t2)} 


(48) 


Introducing  ^  and  -  y  then  yields 


C  =  H  (r  -  \)  {-  -5  [1  -  e'^(T'X)j  +  x  e"^(T'v)  +  O  (r2)} 


(49) 


This  response  is  illustrated  in  Fig.  4  and  may  be  somewhat  more 
realistic  than  the  previous  models.  Again,  the  peak  height  is  C  =  T. 


Reverting  to  physical  variables  where  Pq  is  the  peak  pressure  and  /3 
is  written  as  T/Tq,  where  T  is  the  time  scale  of  the  pressure  decay, 
gives 


P  at 

r  0 

'E  ^  pU 


P  aT 

&  00 

%  %  pU 


(50) 
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4.  THE  CASE  OF  A  CONSTANT- DEPTH  CHANNEL 


The  governing  equation  for  the  disturbance  in  a  constant-depth  channel 
with  a  wall  at  the  leading  end  follows  immediately  from  Eq.  7 


hr)  -  -  TV. 
xx  g  tt 


—  P 
Og  xx 


The  boundary  conditions  may  be  written  as 


x  =  0 


I  .  ^0 

.  u  =  — 

nx 


x  —  »  for  t  <  •» 


while  the  initial  conditions  are 


T)  =  0  =  T)  at 


t  =  0 


For  small  time,  an  approximate  result  is  again  found  by  neglecting 
gravity  forces  so  that  Eq.  5  1  becomes 

T)  =  -  P 
'tt  p  xx 


Writing  the  pressure  as 


p  =  P<t-$> 


gives 


3  =  — —  P 

xx  tt 


Letting  P  be  a  step  function 


P  =  PqH  (t  -  £) 


where  Pq  is  the  magnitude  of  the  pressure  step,  we  have 

Pxx  =  (56) 

where  6  is  the  delta  function.  Then  it  follows  that 

”tt  =  (t  - (S7) 
n  u 

Integrating  with  respect  to  t  from  0  to  t  and  utilizing  the  initial  conditions 
gives 


Tjtx,  t)  -  ,  {h  |t  -  |j)  -  H  (-  L,  y 

pU 


.  H  tfifeH 


On  the  other  hand,  if  we  take  P  to  be  a  smooth  function,  we  have 


1  P 

«  CL2  " 


so  that 


rj(x,  t)  =  {p  (x,  t)  -  P  (x,  0)  -  t  Pt  (x,  0)} 


For  example,  letting  P  be  given  by 


„  U  ,  r  'x  M 
P  =  Po  trie  {a  lv-jj  -  t j; 


where  erfc  denotes  the  complementary  error  function,  yields 


v  (x’ 0  =  {erfc  fe  '  l)j  '  erfc  (a  V 

pu 


2a  t  f  2  /x  |2TI 

•777  expLa  lu/1  jj 

Figure  5  illustrates  such  a  response. 
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5.  RESPONSE  AT  LARGE  TIME 


It  is  evident  that  the  present  model  can  only  poorly  describe  the  large 
time  response  since  at  large  time  the  pressure  front  will  be  acting  in 
water  which  is  no  longer  shallow.  Furthermore,  as  the  blast  wave 
progresses,  its  strength  and  velocity  eventually  decrease  to  those  of 
an  acoustic  wave,  so  that  our  assumptions  are  invalid.  But,  still,  it 
is  of  interest  to  see  what  the  results  are. 


Denoting  the  Laplace  inversion  integral  of  Eq.  34  by  I  and  making  the 
change  of  variable 


a  =  s  ©  0  =  yr 


we  have 


I 


c+i  * 
c-i  00 


e0(A  -k2/4A) 
A2  +  k2  Ql/4 


dA 


This  is  of  the  form 


(63) 


(64) 


e0f(A) 

r 


G  (A,  0)  dA 


for  which  the  method  of  steepest  descent  is  suitable.  When  Q  »  1 
(r  »  1)  one  may  evaluate  the  integral  approximately  by  deforming  the 
contour  in  the  A  plane  to  pass  the  saddle- points  located  at 


A 


ik 
T  ’ 


ik 

T 


The  steepest  path  should  be  directed  at  angles  3  ff/4  and  jr/4  with  the 

?  2  2  - 1 

positive  real  axis  in  the  complex  A  plane.  Since  G  =  (A  +  k  Q  / 4) 
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has  two  simple  poles  at  A  =  +  ik0/2  it  is  evident  that  the  residues 
must  be  accounted  for  in  changing  the  contour  from  T  to  T  1  (Fig.  6). 
Following  standard  formulas,  we  have: 


I  =  .  *  ■  C  =  »  *  .  C  +  Residues 
^  If  i  J  ji  ^  TT  i  j  ji  | 


^saddle  points  *  *poles 


1  fW  k"3/2 

as  tv  —  —  —•  i  (cos  k  0  -  sin  k  0) 
2*7 T  T  -  i 


,1  r  +  1  ,  , 

+  ie  Slnl~ r~  k  r  ^  1 


(65) 


Equation  65  may  be  substituted  in  Eq.  34  to  yield 


.  r1/4  i 

Cd  7T  .'0 


(cos  k0  -  sink0)  J  (k§)  — — ~ 

° 


co 

+  ^  sin  |k  — J  JQ  (k  ?  )  dk 


0 

Both  integrals  in  this  expression  may  be  explicitly  evaluated: 


CO 

(Cd)  =  \  J0  (k5)  dk 

poles  O  r> 


■  "  x>(^) 


(66) 


(67) 
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3 .  2  Sample  Analysis 


The  value  of  a  sample  is  only  as  good  as  the  water  in  it  can  be  meas¬ 
ured.  This  measurement  is  a  function  of  how  the  water  is  removed 
from  the  sample  chamber  walls,  the  instrument  used  for  analysis,  and 
how  accurately  the  instrument  is  calibrated.  Although  present  pro¬ 
cedures  may  be  adequate  and  possibly  even  optimum,  detailed  quan¬ 
titative  investigations  of  many  aspects  of  the  procedures  would  provide 
a  sounder  basis  for  analysis  results  and  may  well  show  how  improvements 
can  be  made. 


3.  2.  1  Water  Removal 


Experiments  on  the  sampler  so  far  ha\  e  shown  that  water  "sticks"  to 
the  sample  chamber  walls  and  that  heating  drives  the  water  off.  The 
more  water  that  sticks  to  the  walls,  the  more  margin  for  error  exists 
in  the  analysis  results  since  even  small  percentage  differences  in  the 
amount  that  sticks  in  the  preflight  evacuation  and  the  amount  that  sticks 
in  the  postflight  analysis  might  mean  large  amounts  of  water.  A  pro¬ 
gram  to  quantitative!  I\  determine  the  effects  of  the  following  variables 
on  water  retention  is  proposed. 

a)  Materials.  Present  materials  plus  any  others  that  a  litera¬ 
ture  survey  indicates  might  be  good. 

b)  Chamber  Temperature 

c)  Chamber  Pressure 

d)  Water  Concentration 

Suitable  chambers  of  candidate  materials  would  be  tested  on  the  pres¬ 
ently  used  mass  spectrometer  (unless  a  better  instrument  is  found,  as 
discussed  in  the  following  section).  Chamber  surface  area  and  volume 
would  be  comparable  to  that  ot  an  actual  sample  chamber.  Tests  would 
be  run  at  at  least  three  values  of  each  of  the  other  three  variables  so  that 


2S 


and 


The  complete  expression  for  £  then  is 

C  =  -y  +  (Cd>  +(Cd>  (69) 

saddle  points  poles 

Because  of  the  discontinuities,  it  is  convenient  to  indicate  the  effective 
region  of  each  component  of  this  solution  in  a  \  -  r  diagram,  shown  in 
Fig.  7. 

It  is  seen  from  Fig.  7  that  two  fronts  exist,  one  moving  with  the  pressure 
front  as  represented  by  the  ray  X  =  t,  and  the  other  moving  with  ever 
increasing  speed  ahead  of  the  pressure  front.  The  variable  speed  of 
the  outer  front  is  given  by 


$  =  T-V-  -  VT  (70a) 

or  in  physical  variables 

^  (70b) 
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which  is  precisely  the  local  wave  speed  ~\J  gh  since  the  depth,  h  is 
ax.  The  wave  height  is  infinite  there  (although  integrably  singular)  and 
the  surface  in  the  neighborhood  is  essentially  given  by  Eq.  67.  Physically, 
the  present  model  would  predict  a  splash-like  wave-front  driven  forward 
by  the  pressure  front  and  traveling  at  the  local  wave  speed. 

The  decreasing  surface  height  as  \  -  »  is  entirely  due  to  the  contribu¬ 
tion  of  the  saddle  points,  Eq.  68.  By  using  the  convergent  expansion 


F  (a,  b,  c,  z) 


n£  1 


<a>n<b>n 

"  TcT 


(71) 


Eq.  68  becomes 


saddle  points 


(r/X) 

Vx(T 


1/4 


■^TT 


where 


(72) 


A 


r (i/4) 


(73) 


and  t/X  is  small.  Near  the  pressure  front,  r/ \  sc  1,  the  following 
expansion  may  be  used: 


F  (a,  b,  a  +  b,  z) 


r  (a  +  b)  (a)ntb)n 

■  r^T  r  (b)  L  ~\2 

n>  0  [n-> 


.  (1  -  z)n{log  (1  -  z)  -  2y(n  +  1) 
+  0  (a  +  n)  +  0  (b  +  n)  ] 


(74) 
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Then,  recalling  that  ij)(l  -  m)  =  0(m)  +  n  cot  17m,  it  is  found  that  at 
the  front,  r/  X  =  1 » 

<Cd>  =  H-r 

saddle  points 

Combining  these  results,  one  arrives  at  the  sketch  shown  in  Fig.  8. 


(75) 
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6.  DISCUSSION 


In  this  analysis  of  a  blast  wave  passing  a  shore  two  peculiarities  arise: 
1)  the  solution  is  not  valid  near  the  shoreline;  2)  the  free  surface  rises 
under  a'  positive  pressure.  These  are  related  and  may  be  discussed  in 
terms  of  a  constant-depth  channel  which  was  treated  in  Section  4. 


First,  it  is  clear  that  the  disturbance  must  occur  wholly  beneath  the 
pressure  distribution  as  long  as  the  channel  depth  is  small  enough  that 
the  wave  velocity  is  less  than  the  pressure  velocity.  Lamb  (1932)  found 
that  in  the  case  of  an  infinitely  long  channel  of  depth  h,  a  distributed 
pressure  traveling  at  speed  U  will  result  in  a  surface  change  in  phase 
with  the  pressure  if  U  <  and  oppositive  in  phase  if  U  >  yih. 

This  may  be  seen  as  follows.  Choosing  a  coordinate  system  moving  with 
the  pressure,  one  obtains  a  steady  state.  From  the  Bernoulli  equation 

^ -  +  h(x)  =  const 


it  is  found  by  differentiation  that 


1  dP 

p  -ar  + 


dh 

H3Z 


0 


dh 

cTx 


1  dP 


Q 

TT 


u 


Hence,  dh/dx  and  dP/dx  are  of  \  samf  Isign  if  U  >  ~\!  gh.  It  is 
then  easy  to  see  that  a  positive  pressure  gives  rise  to  an  elevation  or  a 
depression  according  to  whether  the  flow  is  supercritical  or  subcritical. 
Since  the  flow  velocity  is  just  the  velocity  of  the  pressure  in  a  stationary 
coordinate  system,  we  have  Lamb's  result. 
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In  the  present  case  of  a  senii-inlinite  channel,  part  of  the  surface  is 
elevated  because  U  >  Vgh.  However,  due  to  conservation  of  mass, 
this  fluid  musi  be  supplied  from  behind  the  pressure  front.  This  is 
in  accordance  vith  the  results  presented  earlier.  The  same  con¬ 
siderations  apply  to  the  case  of  a  sloping  beach.  The  immediate 
neighborhood  of  the  shore  deserves  further  study. 
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PART  II 


GROUND  WATER  FABLE  MOTION 


0 


1.  INTRODUCTION 


In  the  event  of  a  nuclear  explosion  near  the  ground  surface,  a  cavity 
would  be  produced.  If  the  ground  water  level  at  the  location  of  the  nuclear 
explosion  is  relatively  high,  this  would  influence  the  flow  of  ground  water. 
In  particular,  radioactive  contamination  of  the  ground  water  basin  is 
possible.  Previous  investigation  into  deep-underground  explosions 
indicates  that  water  contamination  is  of  little  importance.  The  same, 
however,  may  not  be  true  for  a  near- surface  explosion. 

Several  problems  are  postulated  within  this  report.  These  problems 
are  either  completely  solved  or  solved  to  the  point  that  numerical  eval¬ 
uation  becomes  straightforward.  These  problems  are: 

a)  The  filling  of  a  cylindrical  crater  in  an  unconfined  aquifer. 

b)  The  filling  of  a  cylindrical  crater  in  an  unconfined  aquifer 
near  a  river. 

c)  The  discharge  of  a  stream  into  a  large  well  in  a  confined 
aquifer. 

d)  A  quiescent  well  in  a  uniform  flow.  The  advection  and 
diffusion  of  contaminants  initially  present  in  the  well  are 
also  investigated. 

e)  Flow  into  a  water  surrounded  crater. 


PRECEDING 
PAGE  BLANK 
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2.  THE  FILLING  OF  A  CYLINDRICAL  CRATER  IN  AN 
UNCONFINED  AQUIFER 


The  problem  under  consideration  is  the  flow  of  ground  water  into  an 
empty  cylindrical  reservoir.  The  physical  occurrence  of  this  problem 
may  be  the  filling  of  a  cylindrical  void  excavated  by  a  nuclear  explosion. 
This  situation  is  shown  in  Fig.  9. 


GROUND 


JO. 


ER  TABLE 


“7V/T7  ~/~7  //////  77T77 

BEDROCK 

BEFORE  EXPLOSION 


POROUS  MEDIUM 

7J7  ////  7 


Figure  9 


The  material  within  a  radius  R  of  the  explosion  is  assumed  completely 
removed.  The  water  in  the  ground  starts  to  flow  out  of  the  wetted  face 
of  the  walls  of  the  crater  and  begins  to  fill  the  crater.  As  time  approaches 
infinity,  the  crater  would  be  filled  with  water  to  the  level  H.  From  the 
mathematical  point  of  view,  this  problem  is  one  of  unsteady,  unconfined 
flow  radially  towards  a  large  "well"  which  is  initially  empty. 
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The  problem  is  first  formulated  by  writing  down  the  differential  equation 
and  the  boundary  and  initial  conditions. 


Z 


R 


Figure  1  0 


The  differential  equation  is  given  by  combining  the  equation  of  continuity 

H  =  o 

~  r  r  '■/, 

with  Darcy's  law 

u  =  ^o,  O  =  -  k  flp/y)  +  z] 

where  axial  symmetry  (  /  ^ 0  =  0)  has  Keen  assumed  and 

r  =  radial  coordinate 
z  =  vertical  coordinate 

u  =  velocity  component  in  the  r  direction 
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w  =  velocity  component  in  the  z  direction 
u  =  velocity  vector 
< o  =  velocity  potential 
p  =  pressure 
y  =  specific  weight 


This  gives 


=  0 
oz 

This  equation  applies  within  the  flow  field  which  is  bounded  by  the  bottom 
bedrock,  the  top  free  surface,  the  left  seepage  face,  and  the  height  H  at 
r  =  infinity  as  shown  in  Fig.  10. 

The  initial  condition  is  the  quiescent  state  with  the  free  surface  given  by 
z  =  H  -  constant  for  R  r  < 


2 

<7  f 5 


l  >  f 

=  0  t'TXVW 


The  boundary  condition  at  the  right  side  ( r  — •  "*=)  is  the  velocity  w  -•  0. 
The  boundary  condition  at  the  bottom  is  w  =  0  (z  =  0).  The  boundary 
condition  at  the  seepage  face  r  =  R  (0  z  ■'  h  )  is  p  =  0  which  implies 
<p  =  -  kz. 


The  boundary  condition  at  the  free  surface  z  -  h(r,  t)  is 


a)  kinematic 


m 


oh  ^  'Hp  _2_h 
t  "Sr 


^  r  *z 


b)  dynamics 


<p  -  -  kz 

where  m  is  the  porosity. 
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We  now  nondimensionalize  the  variables  by  the  following  formulas.  The 
nondimensionalized  variables  are  shown  in  Fig.  11. 


h*  =  h/H,  z*  -  z/H,  r*  =  r/R, 


<p*  =  g)/kH, 


t* 


Hm7k 


Figure  11 


The  two  boundary  conditions  on  the  free  surface  may  be  combined. 
Ignoring  the  stars,  we  have 

O  -  -h  -  o(  r,  /  h(  r ,  t),  t )  -  -h(r,t) 

Differentiating 

^h  _  {  oo  £|h  ^  oh  (  h(p\ 

^  t  v  t  Az  6  t  -  t  1 1  \  dz J 
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Sh  _  S<p  So  _dh  ^  Sq  _  .Sh  f  ^ 

Sr  Sr  Sz  Sr  Sr  Sr  \  Sz/ 

These  may  be  used  to  eliminate  one  of  the  variables  from  the  S^  (free 
surface)  conditions. 

2.  1  Method  I:  Expansion-in-time  Series 
- - 

! 

I 

One  of  the  standard  ways  of  attacking  this  type  of  problem  (initial  boundary 
problem)  is  the  expansion  in  powers  of  t.  Thus,  we  let 

1 

2 

O  (r,  z,  t)  -  Oq  (r,  z)  +  tOj  (r,  z)  +  t  +  •  •  • 

h  (  r,  t)  =  hQ(  r)  t  th  ^  (  r)  i  ... 

Substituting  into  the  differential  equation 

r  — -  0  i  =  0,1,2 _ 

p  i  H/R  =  constant 
The  boundary  conditions  are 

r  -  30  i  =  0,  1,  2,  .  .  . 

z  =  0  i  =  0,  1 ,  2,  .  .  . 

Qj  =  0  i  =  1,2,...  on  r=l,0<z<l 

(pQ  =  -  1 ,  onz=l  l<r,c® 
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SO. 

i 

Sz 


=  0 


d"V 


oz 


3  0 


<Pn  =  -z 


°i  ,21  Jj_  ' 
2  P  r  or  1 


oz 


For  <p.  i  =  1,  2,  .  .  .  on  S^.,  the  condition  will  have  to  be  derived  by 
substituting  the  series  into  the  boundary  condition.  The  boundary  con¬ 
dition  is 


do 

2  do  7r  do 

P  dr  '  dp  d z 

dz 


0 


or 


do  ,  2  (  doN  ‘ 

-*p  K~ / 


.70 

■§7 


/  -o' L 


=  0 


Substituting  the  series  into  the  boundary  condition, 

-i  2  .. 


o,  +  2«j2+  ...  «  »-  Lo0  .  If,  i 


,  1  toi 

/.  7. 


+  o0  lio1 

A  7, 


Separating  the  various  order  of  t 

°l  *  •  {<>  KW)  *  \TTJ  '  {—>  ■  on  2  "  1 


‘j 

J 


o 


2 


o  r 


\ 

i 


e  r 


dz  d  z  J 


etc . 

The  advantage  of  this  method  is,  as  is  llu*  case  with  most  approximate 
theories,  that  the  nonlinear  system  is  reduced  to  a  succession  of  linear 
problems.  The  solutions  should  be  valid  for  a  relatively  small  t. 
Remembering  that  it  is  the  dimensionless  time  which  must  be  small,  it 
is  instructive  to  observe  that 
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Typical  values  of  H,  m,  and  k  indicate  that,  for  t*  =  0.  1,  t  is  of  the 
order  of  a  day  or  so.  Thus,  the  small  time  fcf  solution  is  not  entirely 
useless.  Indeed,  the  solution  </3q  +  t<Pj  should  be  good  for  a  period  of 
a  week  or  so  after  initial  blasting. 

The  problem  for  </5n(r,  z)  is,  therefore,  as  follows  (Fig.  12); 


4>  -- 1 


*=-z 


t -r  +  p27^(: 


IT*  °-  =  _1 


Figure  12 


Before  attempting  a  solution,  let  Oq  -  -  1  +  Oq  so  that  the  problem  in 
3^  is  (Fig.  13) 


*1=0 


#>I-Z 


2  / 

5  2  i  a  /  . 

3z2  p  7  a r{rTr)  -  0 


<K  s  o 


Figure  13 
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To  obtain  <p q  we  separate  variables. 


<Pq  =  Z(z)  R(r)  to  obtain 


Z '  2  1  d  /  dR\  4  ,2 

T  =  7R  =  constant  =  -k 


Z  *  k  Z  =  0 


=  { 


sin  kz 


cos  kz 


J 


°°0 

-  -  0  at  ;',I  0  :Z  *  A.  coskz 

5  z  k 


=  0  at  z  —  1  I  k (n  I  -)  77,  n  •  0,  1,  2, 


The  equation  in  R  is 


=  0 


d^R  1  dR 

.2  r  dr 

dr 


(k Ip)1  R  -  0 
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which  gives 


R  =  Cl  (-  r)  +  D  K  (-  r) 

op  op 

where  I  and  K  are  the  modified  Bessel  functions  of  the  1st  and  2°^ 
o  o 

kinds,  respectively. 


The  boundary  condition  for  r  —  00  (<Pq  -*  0)  then  requires  C  =  0  since 


and 


K  (x) 
o 


e 


x 


Thus,  solution  may  be  expressed 


O 


0 


n=0 


K 

o 


«" 1  i)v 


It  remains  to  evaluate  which  may  be  obtained  by  the  use  of  the  re¬ 

maining  boundary  condition,  i.e., 

at  r  -  1 ,  Oq  =  -  z 

r‘-  An  +  i 

L  An  Ko  (  ~p  J  COS  (n  +  l)lTz  =  l'z 

n=0 


Let 


r(n  +  2)ffi 


B  =  A  K  .  -  . 

n  n  o  l  p  J 
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Then,  by  Fourier  series  methods,  Bn  is  easily  found  to  be 


8 

nZ(2n  +  1)* 


Hence, 


K 


<Pr 


1  + 


\  «  _ 

^n772(2n+l)2  K  , - j 

n=0  o  V  2 


/  2n  +  1  1Tr\ 

o  V  2  O  /  i  _  z 

/?„  n  <FY  c°sUnt  11  "I 


Thus,  we  now  know  OqU,  z). 


Let 


f(: 


2 

2 

J  *  (  t  g\ 

i  +  -ol 

1  P  \  ^ —  J  1 

1  —  - 

r  \  cr  / 

\  JZ  / 

oz  1 

z-  1 


The  problem  for  is  then  (Fi^.  14) 


Z 


which  may  be  solved  by  Hankel  transform  techniques.  These  will  not  be 
performed  here. 
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2.  2  Method  II:  Quasi-Steady  Solution 


The  method  presented  before  would  yield  the  exact  solution  if  the  series 
were  convergent  and  if  the  entire  series  were  found.  Practically,  this 
should  be  convergent  for  a  small  enough  t  and,  for  such  a  small  t, 
perhaps  one  or  two  terms  in  the  series  would  be  sufficient. 

Even  the  first  two  terms,  however,  are  not  easily  evaluated 
nume rically,  although  it  can  certainly  be  done.  This  numerical  difficulty 
would  be  more  severe  in  the  case  of  than  <Pq.  In  other  words,  even 
though  the  method  of  expansion  in  powers  of  t  yields  analytical  solution 
in  the  form  of  o50(r,  z),  p  j(r,  z),  etc.,  to  obtain  numerical  results,  it  is  still 
necessary  to  do  a  substantial  amount  of  computational  labor.  Therefore, 
it  would  be  desirable  to  have  a  simpler  solution  which  would  give  numerical 
results  without  too  much  work.  This  would,  of  course,  have  to  be  a 
more  approximate  theory.  Such  a  theory  is  presented.  This  theory  is 
based  on  the  assumption  of  quasi- steady  motion  and  the  Dupuit-For scheimer 
approximation. 


Z 


^ _ • _ i 

i 

i 
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D 
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Figure  1  5 


Let 


Q(t)  =  the  discharge  through  the  sides 

h(t)  =  the  height  of  water  in  the  crater 

H  =  the  height  of  water  in  the  porous  medium  for  r 
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— *  CD 


At  a  given  instant  of  time,  it  can  be  assumed  that  the  motion  is  given  by 
the  Duprit-Forscheimer  solution.  With  this  assumption, 

h2  '  n2  =  1FR-  loR  r/R  (1) 

where  K  is  the  permeability,  and  r  is  the  radial  distance  representing 
the  intersection  of  the  Dupuit  free  surface  with  z  =  II. 

There  are  now  three  unknowns  h(t),  Q(t)  and  r(t).  Therefore,  two  more 
equations  are  needed.  One  of  them  is  obviously 
♦ 

!  Q(r)  dT  TTRZh 
*0 


or 


Q 


Tt  R 


>  dli 
d  t 


U-) 


The  second  one  is  not  so  easily  found.  One  way  is  that  the  volume  of 
water  in  tin  reservoir  should  come  out  of  the  voids  of  the  dewatered 
region.  That  is 

l  ■' 

ttR  h  -  e  mp  s (p)  dp  (3) 

'R 

where  s  is  the  drawndown  as  illustrated  in  Fig.  8  and  is  given  by  the 
Dupuit  formula 

s(p)  ,  Vir-^-io^e 

and  f  is  the  porosity. 


48 


Figure  16 


Before  attempting  a  solution,  nondi mensionali  ze  by  the  following  variables 


r*  =  r/R,  h*  =  h/II,  Q  =  Q/*KH2,  t*  =  -J- - 

ir/HK 


Then 

1  -  h*^  I  Q-  log  r* 
dh* 

Q  ::  =  cl  t  * 

,r:;: _ 

h:‘:  =  2e  *  P  (1  -  -\/  1  -  > .  v!'  log  p)'  dp 
“  !  V 

rr-“  1  1  P1*"  i — - — , 

=  —  -  2J  -  **  V  1  -  Q*  log  0  dp 

2  p**  / - - 

-  C  (r*  -  1)  *  I*  J  P  -\/  1  -  Q:|:  log  p  dp 

Dropping  stars, 
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i  j  mini 

'  f  ■  JtnWWwinWWiB 


1  ■  b2  =  log  r 

h  =  €  (r  1}  '  2f  Jj  P\f  "  fr  lo8  P  dp 

This  is  a  system  of  nonlinear  into,-  ro-differential  equations, 
examine  the  integral 


l  - 

1  =  W 1  ■  Q  loK  P  dp 

Let 

P  =  cy  dp  =  t.y  dy 

then 

1  =  j  ^  V" 1  -  Qv  dy 

Now  let 


7. 

1 


1  -  Qy.  y  -  (J  -  Z)l Q,  dy  ft  _  dz/Q 


"“'-o  (1-Z,J  Vs" t  =  z^~  p<4 


Q 


Let 


w 


1 


w  - 


2 

x 


;2/Q  i^y 

V2  / 


5/2 


-  w 


e2/Q 

Q 


x'dx 


First 


d/. 


^0 


r  -x2  2 

The  integral  Je  x  dx  may  be  integrated  by  parts  to  give 


r  -x2  2 . 

e  x  dx 


J  (x/2)  e  X  2x  dx  =  ~  e’X  "  j  j*e  X  dx 


The  original  integral  I  is  evaluated  between  the  limits  r  and  1.  The 
relationship  between  x  and  p  (after  tracing  through  all  the  changes  of 
variables)  is 

x  =  2/Q y"z"  -  -yj  2/C)  "yj  1  -  Qy  =  y2/Q  yi  -  Q  log  p 


Hence 


»r  / _ 

p  yi  -  q  log  p  dp 

iVTvi  •  Q  i°k  r 


x  -x 
2  e 


Cvi  -  Q  log  i 


\^7q 


1  -  Q  log  r  2 

e  dx 


1  -  Q  lot 
2  Q 


r  r  2(1  -Qlogrll  1 

-  e*P  L - q-M  '  Zn 


tJIq 


i  {(  V  2/Q  yi  -  Q  log  r  )  F  -  (  yiTo)  f} 


i**  _  *2 

where  F(x)  =  e  9  d' 

J0 

To  recapitulate,  one  has 


1  -  h  =  Q  log  r 


Q  =  dh/dt 


h  =  C  (r2  -  r  exp  [  -  ~ -  g  '°^  r)] 


■2/Q 


4 


±  •  2g  lo6  r  |  F 


-  (V^Tq)  F 


to  be  solved. 


Now  consider  the  function  in  the  brackets  for  the  expression  in  h. 


F-\f: 


2  -  2Q  log  i 
Q 


-  exp  (-  — 


-  2Q  log  r  ' 


j  -  ’sfllQ  e 


-2/Q 


fy— §i0t!  7  F  4  <V^»F 


Y-2  ■  2g  r  V‘  -Ql°«r 


so  that  the  expression  E  may  be  rewritten 


^Z/Q  ^j  \  -  Q  log  r  exp  ^  (1  -  Q  log  r)J  -  -\jzlQ 

(V2/Q  V1  ■  Q  IoS  r  )  F  4  (-\J2/Q)  F 


-2/Q 


■2/Q 


{  -0  -  Q  log  r  e2  log  r  -  l} 


("\/2/Q  VT-  Q  log  r)  F  +  (YI7Q)  F 


■  V27Q  e-2/Q  {Vi  -  Q  loi!  r  r2  -  l} 
(V27Q"VTTQ1°n)F+  ("\/2/Q)  F 


The  system  of  nonlinear  integro-differential  equations  can  be  solved, 
and  the  three  functions  Q(t),  h{ t)  and  r(t)  may  be  obtained.  This 
involves  some  numerical  computation;  however,  the  amount  oJ  com¬ 
putation  should  be  much  less  than  the  computation  for  <0^  and  Oj,  etc. 
of  the  first  method. 

Another  method  will  be  presented  which  is  not  based  on  Eq.  3  but  is 
based  on  the  formula 


r 


(3a) 


or 


1  +  1 .  5  'yfUlHK  m\fF 


where  v  =  Kb/e.  b  is  a  weighted  average  of  the  depth  of  the  flow 
which  may,  for  all  practical  purposes,  be  taken  as  equal  to  H  . 


Letting  1 . 5  u/HK  =  /3,  we  have  (dropping  stars), 
r  =  1  +  /3 


(3b) 
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Equation  3a  simply  states  that  the  zone  of  influence  of  the  draining  crater 
spreads  at  a  rate  proportional  to  yF  (which  is  normal  for  diffusion- type 
pliehcmena)  with  a  proportionality  constant  equal  to  1. 5  V*  This  latter 
constant  is  based  on  observations  of  discharging  wells*  Hantush  (1965), 


The  other  two  equations,  of  course,  remain  unchanged  and  are 


h2  =  Q  log  r 
Q  =  dh/dt 


repeating  Eq.  3b, 

r  =  i  +  /syr 


These  combine  easily  to  give 


dh/dt  log  (1  t  pyt> 


which  may  be  written 

dh  _  _ dt _ 

1  -  h2  log  (1  4  0^/7) 

*2i. 

The  LHS  is 


which  integrates  to 

'  1  /2 

-  j  lo8  (1  -  h)  t  j  iog  ( 1  +  h)  =  log  (j-7-^) 
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The  RHS  is,  letting  t  =  w  ,  w'  =  /?w,  ?.  =  1  +  w'  successively 


2w  dw  _  2  w '  dw '  _  2  (z  -  1 )  dz 

log  (1  +  log  (1  +  w')  “  log  z 


2  T  /  z  dz\  dz  1 

L  Vlog  v 


L  Uog  zJ  log  z  j 

2  r  d  (z2)  dz  1 
02  L  log  z2  ‘  log  55  1 


f  dx 
J  log  X 


=  ioS|ioS  x|  +  i^i 


Hence,  h(t)  is  given  implicitly  by  the  relation 


log  if^hl  =  lo«l  1°K  U  +  ?  Vt)2|  -  log  |  log  (1  + 

+  log  (1  f  $  ^/l)2  -  log  (1  + 

2 

+  yt~z  { (log  (i  +  £  y^2)  -  r log  ( i  + «  y?)] 


+  ... 


(yr^)  =  exP  (f(x2)  *  fU) 

where  x  =  1  + 


(•X  d! 

J,  loe  ? 


JTT  =  exp  (2  [f(x2)  -  f(x)3 )=  F(x) 
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one  has 


h  =  £j4 
ftt 

The  function  F(x)  and,  hence,  h(x)  may  be  evaluated  simply  by  means 
of  the  digital  computer.  It  should  look  like  the  curve  shown  in  Fig.  9. 


Figiro  17 


* 

The  function  h(x)  is 


h(x) 


eBlxl  .  I 
esUl  1  1 
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where 


x 


This  is  the  same  as 

8<X>  =  -Cob  x  V*  *  Ei  (los  x2)  -  Ei  =0 

Ei  (2  log  x)  -  Ei  (log  x) 
=  Ei Ul)  -  Ei  (l) 


where  "t  =  log  x 


Thus, 


g  ( l )  , 

h(logx)  =  hU)  =  2.  jlA 


+  1 


tanh  r^/Ei(2i).  Ei(^j 


1  =  log  (l  +  yT) 


W)  i. '  -v.lu.Uc 1  in  the  ranse  0  «  t  *  2  and  this  is  shown  Braphica,ly  in 

.8-  10.  I  can  be  seen  that  this  las,  method  ttives  a  (Ulin,  ,ime  which 
depends  only  on  the  Quantity  R  /i/ 

q  amity  k  /y.  The  accuracy  of  the  result  depends 

on  the  formula  r  =  R  +  1. 5  y^tT  P 
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3.  THE  FILLING  OF  A  CRATER  NEAR  A  RIVER  IN  AN  UNCONFINED 


AQUIFER 

The  problem  of  the  flow  of  ground  water  into  an  empty  cylindrical  crater 
in  a  medium  with  infinite  extent  has  been  considered  in  Section  2  and  is 
relatively  simple,  due  to  circular  symmetry.  In  the  case  of  a  crater 
located  near  a  river  or  a  lake,  such  as  Detroit,  the  problem  becomes 
rather  complicated  due  to  its  nonsymmetry.  An  exact  solution  involves 
a  considerable  amount  of  computer  calculation  which  is  beyond  the  time 
limit.  In  order  to  obtain  the  numerical  results,  approximate  solutions 
are  given.  These  approximate  solutions  are  based  on  some  assumptions, 
in  particular  the  Dupuit- Forscheimer  assumption,  as  has  been  discussed 
in  the  previous  sections. 

Figure  19  indicates  the  general  configuration  of  the  problem.  The 
elevation  of  the  water  along  the  river  is  assumed  to  be  constant.  That 
is  to  say  the  fluctuations  of  water  level  of  the  river  are  neglected.  The 
mathematical  formulation  of  the  problem  is  similar  to  the  problem 
indicated  in  Section  2  except  along  the  river  where  p  =  k(p/r  +•  z)  is 
constant.  The  differential  equation  and  boundary  conditions  can  be 
written  as  follows: 

The  differential  equation 

V2<p  =  0 

The  boundary  conditions  at  the  free  surface  are  the  same  as  in  the  case 
of  an  infinite  aquifer,  that  is 

a)  Kinematic 


r 


ELEVATION  VIEW 


b)  Dynamic 

</3  =  -  kz 

The  boundary  condition  along  the  river  is  <p  =  -  kH. 

The  notations  are  as  defined  in  the  previous  sections  except  X  which 
is  the  distance  between  the  crater  and  the  river. 

One  of  the  methods  of  solution  to  this  problem  is  to  use  the  method  of 
images.  On  the  opposite  side  of  the  river,  as  shown  in  Fig.  20,  one 
introduces  a  recharge  crater,  with  the  recharge  strength  equal  to  the 
flow  strength  of  the  crater.  One  obtains  the  solution  for  each  individual 
crater  by  use  of  power  series  as  has  been  done  .‘n  the  previous  section. 

Since  the  equation  of  the  velocity  potential  for  different  order  is  linear, 
we  may  superpose  the  solution  of  two  individual  craters  and  obtain  the 
resultant  solution  which  satisfies  the  boundary  condition  along  the  river*. 
As  has  been  demonstrated  in  the  previous  sections,  the  calculation  of 
the  numerical  solution  of  one  crater  by  use  of  power  series  expansion 
is  quite  laborious.  Therefore,  in  the  following  section,  an  approximate 
solution  is  given  from  which  numerical  results  may  be  calculated  without 
too  much  work. 

It  is  clear  that  the  solution  of  this  problem  can  be  divided  into  two  stages. 
Stage  1,  when  the  radius  of  influence  zone,  r,  is  less  than  the  distance 
between  crater  and  river,  X.  In  this  stage,  the  problem  can  be  consid¬ 
ered  as  if  the  river  were  not  present.  The  solution  has  been  obtained  in 
the  previous  section.  Stage  2,  when  the  radius  of  influence  zone,  r, 
is  greater  than  X.  In  this  stage,  the  river  begins  to  effect  the  motion  of 
the  ground  water.  However,  in  this  stage  the  gradient  of  the  velocity 
potential  tp  becomes  smaller  (the  gradient  of  velocity  potential  is  max¬ 
imum  at  initial  instants  and  near  the  crater).  The  assumption  of  small 
drawdown  may  be  used. 
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As  the  drawdown  becomes  smaller  the  discharge  in  a  column  of  unit 
width,  q,  may  be  approximated  as 

q  =  -  kHvh 


where  H  is  average  depth  of  the  aquifer.  Substituting  this  equation  into 
the  unsteady  continuity  equation 


V  •  q  = 


€ 


3h 

at 


where  €  is  the  norosity  gives 


V2h 


C 

kTT 


3h 
3  t 


This  is  the  usual  heat  equation.  Since  the  equation  is  linear,  one  could 
also  use  the  superposition  method  to  satisfy  the  boundary  condition  along 
the  river. 

Writing  the  Laplace  equation  in  cylindrical  coordinates  one  has 

o2h  ,  J.  oh  _  C  3h 

-  2  r  ar  "  KH"  3t 
or 


The  solution  to  this  equation  with  the  boundary  and  initial  conditions 
prescribed  in  Section  2  is 


h 


_  Q 
°  4TTKH 


-_r_i 

4trHt 


du 


Q 

477KH 
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where 

r2< 

X  =  - — 

4KHt 

The  drawdown  is 


H  -  h  =  - 
°  4ffKH 


—  E.  (-x) 

TT  L  1 


Superposing  two  craters  in  space  as  indicated  in  Fig.  2  l  one  has 


H  -  h  = 
o 


4ffKH 


E. 


2  \ 
rl  C 


4KHt 


+  E. 


2  \ 
r  f 
2  e 


4  K II  t 


1J 


where  r^  and  are  distances  from  crater  and  image  crater, 
They  are  equal  to 


2  ,  »  .2  2 

rj  =  (x  -  A.)  +  y 


2  /  x  2 

r^  =  (x  +  X)  +  y 


in  a  rectangular  system. 


The  equation  governing  the  filling  of  the  crater  is 


Q  =  ^ar 


Substituting  into  the  previous  equation  wc  have 


4HK(H-h) 


dt 

2  \ 

2 

& 

R2 

F 

rl  € 

+  F 

r2  f 

i 

i  4KTTtJ 

1 

(  4KHtJ 

respectively. 
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mmmm 


4.  DISCHARGE  OF  STREAM  INTO  A  LARGE  WELL  IN  A 


CONFINED  AQUIFER 

In  the  case  where  the  explosion  has  created  a  cavity  near  a  stream  and 
the  aquifer  under  conside  1  ation  is  confined,  the  solution  of  the  problem 
of  filling  the  cavity  may  be  rather  simply  found  with  the  notation  shown 
in  Fig,  22. 


Figure  22 


b  =  the  thickness  of  the  aquifer 

Xj  =  distance  from  origin  to  center  of  point  well 

hQ  =  water  level  in  well 

h  =  water  level  in  stream 
s 

d  =  distance  from  center  of  real  well  to  stream 
R  =  radius  of  well 
x  =  coordinate 
y  =  coordinate 
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The  solution  may  be  written  as 


/  X  2 

„  (x  -  x.)  +  y 

h  =  WT  in - 2 - 2 

4  T  (x  +  ^ 


where 


Q  =  discharge 

T  =  hydraulic  transmissivity  =  Kb 
K  =  permeability 


This  is  the  solution  for  a  point  well  at  x  =  Xj,  y  =  0  with  a  stream 
of  head  h  along  x  =  0.  The  distance  d  and  the  well  radius  R  are 

9 

as  yet  unspecified. 
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For  a  finite  diameter  well,  any  equipotential  may  be  taken  as  a  well 

cross  section.  Thus,  if  h  is  the  level  of  water  in  the  well,  then 

w 

along  the  well  edges 

i  v2  2 

q  (x  -  x  )  +  y 

h  -  h  =  tn - U; - 5- 

s  w  4nT  ,  ,  .2,2 

(x  +  xp  +  y 


(x-xi)2  +  ''2  r4rr  „  ,  ,n  . 

- 5 - 5T  =  exp  -pr-  h  -  h  /  .  =  Of 

,  ,  .2.2  _  U  s  w  J 

(x  +  Xj  +  y 


This  is  the  equation  of  a  circle  (the  well  edge), 


(x  -  x^)2  4  y2  =  a  (x  f  Xj  )2  +  ay2 


2  1  +  O  ,  ,2,2 

x  -  7  —  +  Xj  4  y  =  o 


,  (  l  +  oA  ,  (  l  +  o^2  2  f  l  +  oA2  2 

■ 2  (xi  i  -  J  x  vxi  i  -  J  xi  ■  vxi  +  y  -  0 


f  1  +  0/Y  ,  2  2  n  +  ay 

\  1  1  -  (yy  7  1  \  1  -  O'/ 


Thus,  the  center  is  at  x  =  x^  j  ^  ,  and  the  radius  is  Xj  V(h*); 


-  1  =  R. 


This  solution  may  be  easily  adapted  to  a  finite  circular  well  at  a  distance 
d  from  the  stream.  The  center  is  at 


1  +  a  1  +  a 

X1  1  -  O'  *  "  X1  a  -  1 


since  Of  >  1 ,  always,  for  a  flow  into  well. 
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R  =  (tanh  R)  (cosech  ft)  d  = 


=  d  sech  ft  =  R 


The  procedure  for  obtaining  the  solution  is  (Fig.  23): 


y 


Given  R,  d,  h  ,  h  ,  and  K,  b,  etc. 
s  o 

1)  Obtain  ft  from  d  =  (cosh  ft)  R 

d  >  R  ^  =  cosh  ft 

2)  Obtain  from  |x^|  =  d  tanh  ft 

3)  Q  =  2nT/ft  (h  -  h  ) 

s  w 

Solution: 


4ffT 

Q 


h) 


•Cn 


(x  -  d  tanh  ft 
(x  +  d  tanh  ft 


which  gives  equipotential  lines  and,  by  implication,  flow  lines. 
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To  get  filling  times,  etc.  ,  hw  =  hw(t) 


Q(t) 


21TT 


2 ITT  h 


cosh  *  ^ 


rrr  <hs  -  hw<'»  ■  — yh  -  c<tf 


cosh  ■= 

I\ 


c<t)  = 


h  (t) 


w 


One  either  has 


nt 


ffRch  (t) 

w 


Q(r)  dr 


which  gives 


or 


irR  h  C(t)  =  1  Q(t)  cIt 

s 


dC 

dt 


2T 


p  2  .-Id 
R  cosh  -5 
K 


1  -  C  =  c  (1  -  0 


C  =  constant  = 


2T 


p  2  .-Id 

R  cosli 


d  C 

Trz 


=  C  dt 


Integrating 


-  log  (1  -  O  -  C  t  +  constant 
The  condition  t  =  0,  C  -  0  implies  that  the  constant  =  0 


1  -  C 


exp 


2T  t 

R^  cosh  * 
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or 


and  finally, 


h  (t) 
w 


UaP(-^T 3-) 

R 


hw(t)  .  hs[l-exp(-  *11 


R  cosh' 


rr 
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5.  A  QUIESCENT  WELL  IN  A  UNIFORM  FLOW 


If,  after  the  crater  is  filled,  one  can  assume  there  is  a  uniform  flow  over 
the  well  as  illustrated  in  Fie;.  24,  some  fluid  will  pass  through  the  well, 
and  some  fluid  will  bypass  it.  If  the  fluid  in  the  well  were  radioactive, 
and  if  it  were  assumed  that  the  radioactivity  is  simply  advected,  it  would 
be  of  interest  to  know  the  width  of  the  zone  of  contamination  downstream. 


Figure  24 


It  will  be  shown  in  this  section  that  d  =  4  R. 


With  the  x,  y  coordinate  system  shown  in  Fig.  24  ,  the  mathematical 
problem  is 


V  h  =  0  for  r 


+  y  ^  >  R 


-  k  =  U  for  r  -  05 
ox 


h  =  constant  on  r  =  R 
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Without  loss  of  generality,  this  constant  can  be  taken  to  be  zero.  The 
solution  to  this  problem  can  be  found  easily,  and  it  is 


Thus,  the  critical  streamline  (Fig.  2^)  is  represented  bv  the  equation 


as  x  -  r ,  this  implies  \  2R  d/2.  Thus,  the  thickness  of  the  zone  of 
contamination  d  is  4R. 
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Figure  25 


The  foregoing  analysis  is  based  on  nondispersive  flow.  Actually,  the 
radioactivity  would,  of  course,  disperse  transversely  from  the  zone  of 
contamination  with  the  result  of  a  larger  zone  than  4R. 

To  investigate  dispersion  in  a  porous  medium,  the  situation  could  be 
idealized  and  the  case  of  a  simple  mixing  zone  (Fig.  26)  is  investigated. 


The  only  variable  is  the  concentration  C,  and  the  governing  equation 
is  the  dispersion  equation 


77 


ac 


ac  a 


+  VJ7  ~ 


(Dxtt)  +  37  (Dyjf) 


which,  in  the  present  case,  since  v  =  0  and  U  =  u,  becomes 


U 


ac  _  d 

d  x  ”  a  x 


[Dxll  j  +  Ty  [Dy37-J 


Here  and  are  the  dispersion  coefficients  in  the  longitudinal  and 
transverse  directions. 

Neglecting  the  term  compared  with  the  term  on  the  ground  of  a 
boundary  layer  (BL)  assumption,  we  have 

is-  •£{'%) 

CX  ey  \  y  ✓ 

where  €  =  D^./U.  Since  U  is  constant,  would  be  expected  to 

be  constant  also.  Hence,  6  is  constant.  Absorb  C  in  y  to  obtain 


1C 

"ex’ 


> 

r  c 


SY 


BC 


x  &  0 


y  -  x.  C  -  0 

y  -  C  =  1 
y  ''  *■ .  -  A  C  /  e  y  0 


For  a  similarity  solution,  let 

r )  ( y /  V*) 

and  c  =  f(  T))\  then 


:  C 

a  x 


1  ~Tf2 

X 


=  f 
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Thus,  DE  becomes 


i  T?  f*  +  f"  =  0 
and  BC  becomes 


f(«)  =  0  f(-»)  =  1 

The  solution  is  C  =  -  l/Z-^W 

33  2 

i  r  - y 

f(t,)  =  +  -  1  e  d* 

V2 


This  is  the  standard  diffusion  type  solution  and  may  be  found  in  textbooks. 
The  value  of  D  needs  to  be  estimated.  In  Ref.  2  a  graph  shows  D  /D 

y  y 

as  a  function  of  A  V/D  where 
P 

A  =  particle  radius  of  the  bed  material 

P 

V  =  velocity  (U  in  this  case) 

D  =  diffusion  coefficient  (molecular) 

For  large  velocities  V,  (or  for  A  V/D~  10)  D  / D  is  about  5  to  10. 

P  y 

Hence,  can  be  taken  to  be  about  D  to  10D  for  practical  purposes. 

The  similarity  variable  is 

rjz  --  | - X— 

yrvTnu 


X. 


2  VDyX/  U 
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Let  y^Q  be  the  y  where  C  =  1/10 


This  means 


f  =  1/10  >  - -  -  1.64 

2  V<Dy/U)  * 

Thus 

V90  3 

The  mixing  zone,  therefore,  grows  like  ~\J  x  with 


#  3  VlOD  /U 

y 


proportionality  constant 


Applied  to  the  present  problem,  the  flow  situation  would  appear  as  shown 
in  Fig .  27 . 


Since  the  substance  undergoing  dispersion  is  radioactive  in  nature,  there 
is  a  natural  rate  of  decay.  Thus,  the  equation  for  dispersion  should  be 
modified  by  a  decay  term.  For  those  isotopes  that  decay  very  slowly  in 
comparison  with  the  motion  in  the  porous  medium,  the  above  nondecay 
analysis  should  apply.  For  those  isotopes  which  decay  reasonably  rapidly, 
the  analysis  of  the  problem  should  be  as  shown  in  Fig,  28. 

Still  considering  the  simple  case  of  a  mixing  zone, 


Figure  28 


the  governing  equation  is 

—  t  kC  -  D  ^ 

:x  v  l 

■  'jy 

where  it  is  assumed  that  the  radioactive  decay  rate  is  proportional  to 
the  concentration  with  k  as  tin-  proportionality  constant.  Dividing  by  U, 


o C 

IX 


The  special  case  C  -  0  gives  the  solution 


C 


C  exP  Tt  x> 


where  Cq  is  the  concentration  at  x  =  0.  This  may  be  applied  to  the 
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core  of  the  zone  of  contamination  in  the  previous  problem.  This  may  now 
be  taken  as  the  boundary  condition  at  y  -♦  -®, 

Nondimensionalizing  the  variables  by 


C*  =  C/Co,  x*  =  (k/U)  X,  y*  -  V(k/D  )  y 


the  following  equation  is  obtained 


|C*  ,  _  d2C=: 

—  I  U  - - T 


with  the  boundary  conditions 


C*  ■*  0  as  y*  -  »,  for  x*  >  0 


C:’:  -  e  X  as  v* 


x*  •  0 


C*  -0  as  \ 


C*  -  e'X"  f 


jC  ■  -x-  -X  ' 

— x  -  -e  f  i  i* 
j  x 


Zr,. 

d  C  _  -x-  r 

- 7  e  f 

dy ^  ' 


the  equation  becomes 


-  f  +  f 


or 


af  a2f 


The  boundary  conditions  become 


f  -•  0  as  y*  -•  30 

}  x  >  0 

f  -  1  as  y*  -  -  x 


Now  note  that  this  is  the  same  problem  as  the  case  of  nondecay  disper¬ 
sion.  Hence,  the  solution  is 


f 


cl* 


Hence , 


C 


C 

o 


exp  (- 


k 

U 


x) 
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6.  FLOW  INTO  A  WATER  SURROUNDED  CRATER 


In  the  case  of  the  New  Orleans  burst,  flooding  could  extend  to  areas  of 
considerable  size,  and  the  water  will  certainly  surround  the  crater  lip. 
Therefore,  the  time  of  filling  of  the  crater  would  be  different  from  the 
case  of  infinite  media. 


If  the  Dupuit- Forchheirner  assumption  is  used  and  initial  transient 
phenomena  are  neglected,  we  immediately  have  the  equation  for  ground 
water  motion 


H 


2 


r 

o 


R 


The  notations  are  indicated  in  Fig.  2'1,  except  Q  and  K  which  are 
defined  in  the  previous  sections  to  be  discharge  and  permeability. 


Another  equation  which  go\erns  the  filling  of  the  crater  is 


r  > 

1  Q(T)  dr  =  77R“h 

“0 


or 


Q 


77R 


dh 

dt 


Combining  the  equation  with  the  equation  of  motion,  we  have 


7TR 


2 


dh 
d  t 


(  ?  2\ 

[H  -  h  J  rrk 

loS  X 


r 


1 


,<]( 

\i 
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PERVIOUS  LAYER 


£-  IMPERVIOUS  BED  h 


Figure  2? 


or 


dh 


K  dt 


u2  ,2  _  r 

H  -  h  „  2  .  o 

R  log  -g- 


After  integrating,  we  obtain 


1  .  ,-lh 

H tanh  71  = 


Kt 


02  ,  o 

R  logx 


The  initial  conditon,  h  =  0  when  t  =  C 


Therefore,  one  has 


h  =  H  tanh 


H  Kt 


„2  i  o 

R  10RT 


r« 


C 


,  implies  c  =  0. 
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